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The method of McCurdy, Baertschy, and Rescigno, J. Phys. B, 37, R137 (2004) is generalized to 
obtain a straightforward, surprisingly accurate, and scalable numerical representation for calculating 
the electronic wave functions of molecules. It uses a basis set of product sine functions arrayed on a 
Cartesian grid, and yields 1 keal/mol precision for valence transition energies with a grid resolution of 
approximately 0.1 bohr. The Coulomb matrix elements are replaced with matrix elements obtained 
from the kinetic energy operator. A resolution-of-the-identity approximation renders the primitive 
one- and two-electron matrix elements diagonal; in other words, the Coulomb operator is local 
with respect to the grid indices. The calculation of contracted two-electron matrix elements among 
orbitals requires only 0(Alog(A)) multiplication operations, not 0(A 4 ), where N is the number of 
basis functions; N = n 3 on cubic grids. The representation not only is numerically expedient, but 
also produces energies and properties superior to those calculated variationally. Absolute energies, 
absorption cross sections, transition energies, and ionization potentials are reported for one- (He + , 
H^), two- (H 2 , He), ten- (CH 4 ) and 56-electron (CsHs) systems. 

PACS numbers: 


I. INTRODUCTION 

The inherent problem in scaling electronic structure 
methods to larger systems is the prohibitive cost of 
storing and transforming two-electron matrix elements, 
which we denote in chemists’ notation 

[ij\kl\ = J J d 3 fid 3 f 2 Xi(A)Xi(A)|-r^-r|Xfc(r2)Xi(r2) 

(1) 

for a basis {xt}- The set of two-electron matrix ele¬ 
ments is a fourth-rank tensor, such that transformations 
of the set require 0(N 4 ) multiplication operations; so¬ 
phisticated methods such as coupled cluster must cope 
with even poorer scaling, 0(N 6 ). There has been much 
work to circumvent this basic problem HHH, especially 
by Martinez and coworkers. 

We describe a basis set method for electronic structure 
motivated by the discrete variable representation (DVR) 
[fHH that is an adaptation of the method in Ref. m to 
Cartesian coordinates. Starting with an evenly spaced 
grid basis means that the number of basis functions is 
large, but also that the rank of the two-electron matrix el¬ 
ement tensor is automatically reduced from four to three 
due to redundancy. Going further, using the generaliza¬ 
tion of Ref. E2 , we obtain a diagonal two-electron matrix 
element tensor - in other words, we further reduce the 
tensor to the minimum rank one, 

[ij\kl\ ~ dijSki Vi- k . (2) 

(In this equation the indices i,j,k , and l each would run 
from 1 to N = n 3 on a cubic grid.) 


Using this resolution-of-the-identity approximation for 
the treatment of the Coulomb potential within the dis¬ 
crete variable representation, and employing established 
Fourier methods for triple Toeplitz linear algebra mm, 
the computation of two-electron matrix elements among 
molecular orbitals takes 0(Nlog(N)) time, not 0(./V 4 ). 
The method is therefore not quite “linear-scaling”, but 
it is numerically exact; it does not involve any truncated 
sums in a multipole expansion, for instance. 

Gaussian basis sets have traditionally been the pre¬ 
ferred single-electron representation for real-space elec¬ 
tronic structure calculations, due to the localized nature 
of these functions and the speed with which matrix ele¬ 
ments among them may be evaluated. Although Gaus- 
sians have been widely successful, they have inherit lim¬ 
itations in their flexibility; in particular, they are unable 
to represent electrons in the continuum, which is neces¬ 
sary for ionization and electron scattering applications. 
Furthermore, it is not always clear exactly how to obtain 
rigorous error bounds of basis set truncation. 

There has recently been an increased interest among 
researchers in the held to develop grid-based methods us¬ 
ing strictly numerical techniques that can handle a wider 
variety of problems and can be subjected to systematic 
error analysis. A thorough review of grid methods in 
electronic structure can be found in [2(| . Some examples 
of grid-based techniques currently in use are finite differ¬ 
ences frB [2j| , finite elements [2J, [25| , and wavelets [26j . 
These methods make the treatment of arbitrary bound¬ 
ary conditions considerably easier than basis set meth¬ 
ods. Another advantage of grid methods is the flexibility 
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allowed in performing calculations on complicated spa¬ 
tial domains. Finite difference methods are limited in 
this regard since they require strictly rectangular meshes, 
whereas finite element methods offer complete freedom in 
choosing a computational mesh. 

Similar to finite element methods, discrete variable 
representation (DVR) methods have characteristics of 
both a basis set method and a grid method in the sense 
that each basis function is localized around a specific 
grid point, and potential functions are evaluated as local 
multiplicative operators on the grid. Many DVR bases 
have appeared in the literature, including those based 
on Bessel functions fl5 l Lagrange polynomials [27l. [28j. 
and sine functions I29M31I ]. as well as multidimensional 
bases 3 and others described in Refs. mu. 

One issue in evaluating potential energy matrix ele¬ 
ments for molecular systems is how to resolve the sin¬ 
gularities that occur in the Coulomb potential terms. A 
number of methods for doing so have appeared in the 
literature, including the use of energy cut-off functions 
[32 , HH and multipole expansions via Legendre poly¬ 
nomials in spherical coordinates [3, HI]- The singu¬ 
lar Coulomb potential cannot be used straightforwardly 
within the DVR approximation, because doing so would 
entail the use of infinite diagonal matrix elements. 

To address this issue, the present method makes use of 
the fact that the Green’s function for the Laplace oper¬ 
ator is the Coulomb potential. In doing so it follows the 
derivation used in Ref. E3- That work presented a treat¬ 
ment in spherical polar coordinates, using the partial- 
wave expansion of the Green’s function, arriving at ex¬ 
pressions for the two electron matrix elements diagonal in 
the radial index, corresponding to an expansion in Gauss- 
Lobatto DVR for the radial degree of freedom. Here we 
do not use the partial-wave expansion and instead treat 
Poisson’s equation in Cartesian coordinates. 

We present the method in the next section, then re¬ 
sults, and finally in the conclusion we speculate about 
possible elaborations to the method that could make it 
more versatile for excited state and time-dependent prob¬ 
lems, and perhaps even competitive with Gaussian basis 
sets for the computation of results requiring chemical ac¬ 
curacy |35 j]. 


The sine function is defined as 


sine 1 ^ . 

1 1 if x = 0 


and an orthonormal basis in one dimension is 


€i( x ) = y=sinc 


X — Xj 


with A the uniform grid spacing, Xi+i = X{ + A. 


(3) 

(4) 


B. Kinetic energy matrix elements 


The kinetic energy matrix elements among these func¬ 
tions are 


tij — 


i£_ 

2 dx 2 

7t 2 /(6A 2 ) if i = j 

(—l) i_J 7 (A 2 (* — j) 2 ) if i^j 


(5) 


Notice that these matrix elements only depend on i — j, 
i.e., t is constant along diagonals, i.e., t is Toeplitz. A 
derivation of these elements is given in Ref. [3- We 
make a three dimensional product basis in the straight¬ 
forward way, 


Xi(x, y,z) = £i 1 (ar )&2 (y)&3 (z) ■ (6) 

The three dimensional kinetic energy is, as usual, 

Tij = Ul,jlSi2,j2Si3j3 + Siljlti2,j2^i3,j3 + ^il,jl^i2,j2ti3,j3 

(V 

Since t and the identity matrix are Toeplitz, T is triple 
Toeplitz, as are the matrix elements of any translation- 
ally invariant operator. We only use explicit vector-index 
notation in Eqs. [G] and [7] and in sections III FI and III Gl 
In the rest of this paper, we use contracted indices, such 
that for a three dimensional basis function yy(r), or ma¬ 
trix element T)j, the index i (or j) represents a single 
integer that runs from 1 to N = n 3 on a cubic grid. 


II. METHOD 
A. Sine Basis 

Sine functions have been used extensively in several ar¬ 
eas of applied mathematics including numerical solutions 
to ordinary and partial differential equations, interpola¬ 
tion and Fourier analysis [36] but were first introduced in 
the context of a DVR basis for sol ving the Schrodinger 
equation by Colbert and Miller in a good descrip¬ 
tion of the sine DVR can also be found in Ref. i37j . Sine 
basis functions have been used in electronic structure in 
Refs. [IMII- 


C. Discrete Variable Representation resolution of 
the identity for two-electron matrix elements 

In the generalization of Ref. na to Cartesian coordi¬ 
nates, there are several simplifications that result from 
the use of sine basis functions. The present method 
and that of Ref. E3 are founded on the replacement 
of Coulomb matrix elements by matrix elements ob¬ 
tained from the kinetic energy operator via a resolution- 
of-the identity approximation invoking Poisson’s equa¬ 
tion. However, for the two-electron matrix elements, it 
is not necessary to introduce the kinetic energy operator 
into the derivation, if the sine DVR is used. Therefore, 
in this section, we provide the simplest derivation of the 
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two-electron matrix elements used in this method, before 
introducing the kinetic energy operator in the sections 
below. 

The method of Ref. d2 uses the fact that the Coulomb 
potential is the Green’s function of the Laplace operator 
to avoid the inherent problem with using the discrete 
variable representation (DVR) approximation for singu¬ 
lar potentials. It results in an expression, Eg. (1X11) . which 
for the sine DVR basis is equivalent to the resolution of 
the identity described in this section: 

[ij\kl\ = 2TrS ij S k i(w i w j )~^T^ 1 

where the w are the quadrature weights - presently, for 
the 3D Cartesian sine basis uniformly equal to A 3 - and 
T _1 is the limit of the inverse of the kinetic energy matrix 
as the size of the basis is taken to be infinity. Because 
the sine basis is complete in momentum space up to a 
cutoff, the matrix element of the matrix inverse is equal 
to the matrix element of the operator inverse, 

(T _1 )ife = [ d 3 r 1 d 3 r 2 X»(n) i 1 - \ Xk(r 2 ) . (8) 

J 27r|ri 2 | 

Because Eq. [8] holds for the sine DVR basis, there is 
no need to introduce the kinetic energy matrix into the 
derivation of the two-electron matrix elements. Our fi¬ 
nal expression for them, Eq. 1121 results simply from a 
resolution-of-the-identity approximation. 

The resolution of the identity makes straightforward 
use of the interpolating property of discrete variable rep¬ 
resentation (DVR) basis functions: each basis function 
belongs to a grid index, and is zero at all the grid points 
other than that corresponding to its own index. This 
“discrete orthogonality” condition j38| is the defining 
property of a DVR basis set. An arbitrary function can 
be expanded easily in such an interpolating basis, 

f(x) « (uh)~ 1/2 Xi(x)f(xi) (9) 

i 

where Wi is the quadrature weight at point <; presently 
the weights are all equal to A for the one-dimensional 
sine DVR and A 3 for the 3D product basis. This may be 
written as a resolution of the identity, 

f( x )~^2\Xi){Xi\f{x) (10) 

i 

where the integral is performed using the underlying 
quadrature, giving 

Xk(r)Xi{r) ~ hi{wk)~ 1,2 Xk{r f ) (H) 

such that the density (a sum of squares of localized ba¬ 
sis functions) is re-expanded as a sum of localized basis 
functions, without the square. The fact that the auxil¬ 
iary basis, the one in which the density is expanded, is 
the same as the basis in which the wave function is re¬ 
solved means that, at least aesthetically, it is the ideal 
resolution of the identity. 


The expression for the two-electron integral is obtained 
simply from Eqs. Q] and [XT] 

[ij\\kl] « A~ 3 Sij6 k i [ d 3 rid 3 r 2 Xi(n)r^—rXk{r 2 ) ( 12 ) 
J \r 12 | 

D. Application of the method of Ref. Hi to 
arbitrary three dimensional discrete variable 
representations 

Here we provide a complete derivation of both the one- 
and two-electron matrix elements that follows the deriva¬ 
tion in Ref. da closely. This method is founded upon the 
observation that the Coulomb potential is the Green’s 
function of the kinetic energy (Laplace) operator, and 
replaces Coulomb matrix elements in a basis with matrix 
elements obtained from the kinetic energy operator in 
the same basis. There are only two significant differences 
between the derivation in Ref. da and this one: one, we 
use the full three-dimensional Green’s function for the 
Laplacian operator, not its partial wave expansion; and 
two, we eliminate the need for an explicit boundary con¬ 
dition. For sine basis functions, the derivation of the 
two-electron matrix elements may be simplified as in the 
section above, but the derivation here is applicable to 
general discrete variable representations in three dimen¬ 
sions and includes both the one- and two-electron matrix 
elements. 

We define 

V kl (n) = [ d 3 f 2 Xk{r 2 )xi(r 2 ) 1 _ (13) 

J |ri-r 2 | 

so that, with reference to Eq. © , we can write 

[ij\kl\ = y Xi(ri)Xj(ri)y kl (ri) • (14) 


Applying the Laplacian to both sides of Ea. (fl3l) results 
in the Poisson equation 


^\y kl (n) =-Airxk{ri)xi{n) ■ 

(15) 

Here we have used the fact that 


G(r,r) - 

47r|r — r \ 

(16) 

is the free space Green’s function for the Laplace operator 
satisfying 

VpG(f, f) = 5(r — r 1 ) 

(17) 


with the boundary condition G(r, r') —> 0 as |f| —> oo. 
We now approximate y kl {fi) as a linear combination of 
the basis functions, 

y kl {ri) ~ ]T y k n l X n(n) , (18) 

ALL n 
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with an infinite sum over the product basis functions cov¬ 
ering all space. In numerical calculations a finite basis 
is always used, but by including all possible product ba¬ 
sis functions in the above sum we avoid the need for a 
boundary condition term, as was needed (and easily han¬ 
dled) in the prior treatment in spherical polar coordinates 
E2- Applying the Laplacian with respect to r\ to this 
expansion gives 

V^ fei (ri)« Y Vn^kXnir i). (19) 

ALL n 

Multiplying Eas. dlSl) and m through by Xm(fi), inte¬ 
grating over r\ and equating the results leads to 

Y yn T mn = 2TT f d 3 rXk(r)xi(r)Xm(r) ( 20 ) 

ALL n 

where T mn = — \{Xm\^ 2 \Xn) are the kinetic energy ma¬ 
trix elements. 

The integral on the right hand side of Eq. (l20l) is evalu¬ 
ated by a resolution of the identity, Eq. [TTJ such that the 
density (a sum of squares of localized basis functions) is 
re-expanded as a sum of localized basis functions, with¬ 
out the square. From Eas. (l20l) and (fill) . 

Y VnTmn = ZTrSklSlmiwkWl)- 1 /' 2 (21) 

ALL n 

giving 

V% =2n5 kl {w k )- 1 {T- 1 ) nk , ( 22 ) 

where Tis the limit of the matrix inverse as the size of 
the matrix goes to infinity. Using these coefficients and 
inserting Ea. (fl8l) into Ea. lfTdll gives 

[ij\M\ = Y y n [ d3f i Xi(ri)xj{n)Xn(n) ■ ( 23 ) 

Once again, the integral on the right hand side is eval¬ 
uated by resolving the identity and employing DVR 
quadrature to obtain the final expression for the two elec¬ 
tron matrix elements 

[ij\kl\ = 2ir6ijdki(w i Wk)~ 1/2 Tff 1 . ( 24 ) 

E. One-electron matrix element 

The expression for the one-electron matrix element for 
a nucleus at position R 1 here denoted [/,;,■, 


within the present method follows simply from the analog 
of Eq.QS 

Uij(R) ~ Y U nXn(R) ; (26) 

ALL n 


Eq. ED which is the resolution of the identity approxi¬ 
mation; and Poisson’s equation: 

S7\U i3 {R) = 4t T X i{R)Xj{R) ~ 4nA- 3 / 2 6 ijXi (R) (27) 

within the weak variational formulation, i.e., multiplying 
by the left by J d 3 R Xk(R) for all k. The one-electron ma¬ 
trix element is thereby simply related to the two-electron 
matrix clement, 

Uij (R) ~ (R) ; (28) 

both are defined via Eqs. [IS] and [22] in terms of the ki¬ 
netic energy matrix elements. However, we find that ac¬ 
curate results are only obtained when the nuclei are at 
positions R coinciding with electronic DVR grid points 
r). In other words, for the moment, the method requires 
that the nuclei be placed on the Cartesian grid points, 
which is a major limitation. We do not understand this 
behavior and comment upon it, and ways to avoid it, in 
the conclusion. 

In order to use the above expressions in a practical way, 
for general DVR bases, we must have a method of ap¬ 
proximating the matrix elements of T _1 . However, sim¬ 
ply inverting the finite-dimensional kinetic energy matrix 
is not a good approximation and would also require the 
storage of a dense matrix, which is impractical for even 
modestly sized grids. We use the method below, which 
directly calculates the entire set of diagonal two-electron 
matrix elements at one time, and which would be ap¬ 
plicable to other generalized DVR basis sets besides the 
Cartesian product sine functions used here. 

F. Kinetic Energy Inverse 

The key to a practical implementation is to consider 
the representation of T _1 on an infinite grid, and then 
truncate the matrix elements to a finite grid. The full 
kinetic energy in three dimensions is given in Eq. [7] Since 
t is Toeplitz, T is Toeplitz with respect to each of the 
spatial indices, i.e. T is triple Toeplitz, and therefore it 
could be denoted using a symbol with only one three- 
vector index, e.g. TV = u^_j. (In this subsection we 
momentarily revert to explicit vector index notation.) Its 
inverse inherits this property, 

(T~\ = v;_j (29) 

The expression TT _1 = 1 can be rewritten 
00 , 00,00 

(u*v)j= Y U i-o V i = ^ 1 , 0 ^ 2 , 0 ^ 3,0 ( 30 ) 

2 —— 00 , 00,00 

where * represents the discrete convolution product. 

For the sine DVR, we may derive an exact expression 
for the matrix element amenable to quadrature, but in¬ 
stead we solve for all of the matrix elements simultane¬ 
ously using the following method, which is also applicable 
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to other bases. The strategy is to take the lowest-order 
Taylor series expression for the matrix element, for those 
matrix elements at long range, and solve for the remain¬ 
der. 

Thus we approximate the matrix elements of T -1 be¬ 
tween basis functions far separated as 

|i|->oo A 3 1 _ A 2 1 
V ? 2^-27^+^ 1 J 


such that [ij\kl] —>■ SijSu^r-^- We assume that Eq. (fTTTT) 

holds exactly for three-indices i in which ij , * 2 , or is 
greater than n sma u , where n sma u is an adjustable param¬ 
eter, and solve Eq. (Idlil) for the remainder of v. In other 
words, we solve 


f^big i^big i^big 

Y u i-] v i = fiji, oSj2, oSj3,o (32) 

^— Ubig i^big l^big 

given = T A and 


A 2 1 
2tt + i\ + j 2 


( \h\ > n Bm aii or \ 

|^2| 'Tlsmall Or I (33) 

| ^3 | ^ Ismail J 


for the remaining (2 n sma u + l) 3 elements of v. The infi¬ 
nite sum in Eq. [30] is truncated at ribig and we therefore 
have two convergence parameters, n sma u and ribig defin¬ 
ing the approximated T~ l . We have chosen a default of 
40 and 240 for these numbers, respectively, and we verify 
the convergence as a function of these parameters of all 
the results presented below. 


G. Triple Toeplitz linear algebra with Fourier 
transforms 

We continue with vector index notation in this section, 
after which we revert to condensed index notation. 

To construct a two-electron matrix element among con¬ 
tracted basis functions 

Mr) = y, c i a Xi(r) (34) 

i 


we must perform the sum 

M7<5] = X! 27r ( w; i w; fc)" 1/2 ^-fe c | a c i / 3 C £ 7 C fc5 ( 35 ) 

ik 

wherein a triple Toeplitz matrix-vector multiplication is 
performed by the triple Toeplitz matrix v upon the den¬ 
sity <j)%4>p to produce a potential that is then integrated 
over the density or vice versa. 

The matrix T~ 1 is triple Toeplitz (a.k.a., 3-level 
Toeplitz), i.e., (r _1 ) iife ,i'j'fc' = fj-i' ,j-j',k-k' where v is 
a (21 — 1) x (2m — 1) x (2n — 1) tensor and T _1 is an 


N x N matrix with N = Imn. A triple Toeplitz matrix 
is (a) block Toeplitz, e.g., 


T~ i 


(T- 1 ) o (T- 1 )! ••• (T -1 ) n 

(T- 1 ).! (T- 1 ) o ••• (T -1 ) n _i 

: (T- 1 )-! i 

(T _i )_ n _! ; ••• (r- 1 )! 

(T -1 )_ n (T _1 )_ n _! ••• (T- 1 )o 


Furthermore, (b) the blocks are double Toeplitz (or 2- 
level Toeplitz), i.e., they are block Toeplitz with Toeplitz 
blocks (also called BTTB in the literature). Similarly, 
a triple circulant matrix C is such that C%jk,i>j'k' = 

C-i—i' (mod (mod m),k—k' (mod n) (^ ^ X TTl X 77, 

tensor, C is an N x N matrix, N = Imn). 

Triple circulant matrices are diagonalized by the three- 
dimensional Fourier transform (Theorem 5.8.4 in 0): 

C = F’ ■ diag(Fc) • F 


An N x N triple Toeplitz matrix can be embedded into 
an 8 N x 8 N triple circulant matrix 0- Therefore, just 
as with single Toeplitz [39}, a matrix-vector product in¬ 
volving a triple Toeplitz matrix, such as that required 
to compute two-electron matrix elements, may be com¬ 
puted in 0(N log N) floating point operations using a 
fast Fourier transform, instead of 0(N 2 ). Memory use 
is minimal, due to the redundancy inherent in Toeplitz 
matrices. 

The embedding that is used 0 to transform the NxN 
triple Toeplitz matrix into an 8 N x 8 N triple circulant 
matrix consists in padding the tensor v with zeros. We 
summarize the algorithm: 


Matrix-vector product y = (T 1 )* with 

T- 1 : N x N Triple Toeplitz matrix, N = Imn, defined by 
v : (21 — 1) X (2m — 1) x (2n — 1) tensor per Eq. [29] 

1 : *2 = 0 (21 x 2 m x 2 n) 

2 : * 2(1 : l, 1 : m, 1 : n) = x(\, :,:) 

3: V 2 = 0 (21 x 2 m x 2 n) 

4: V2(2 : 21,2 : 2m, 2 : 2 n) = v (\,:,:) 

5: f x = FFT-3D(:E2) 

6 : f v = FFT-3D(u 2 ) 

7: f v = f x x f v (element-wise product) 

8 : 2/2 = IFFT-3D(/ y ) 

9: y(:,:,:) = 1 / 2 (l + 1 : 21, m + 1 : 2m, n + 1 : 2 n) 


III. RESULTS 


A. DVR Method vs. Variational Method 

Remarkably, the treatment we have outlined appears 
to perform better than the variational method (in which 
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Energy 

Virial theorem 

State 

DVR 

Variational 

Exact 

DVR 

Var. 

Ex. 

Is 

-1.9765 

-1.9526 

-2 

-0.4939 

-0.4817 

-0.5 

2 P 

-0.4998 

-0.4953 

-0.5 

-0.4998 

-0.5184 

-0.5 

2s 

-0.4976 

-0.4826 

-0.5 

-0.4987 

-0.5280 

-0.5 

CO 

-0.2189 

-0.1761 

-0.22... 

-0.5282 

-0.7108 

-0.5 

3d* 

-0.2155 

-0.1712 

-0.22... 

-0.5494 

-0.7184 

-0.5 


TABLE I: The n = 1, 2, 3 eigenvalues and virial theorem ra¬ 
tios < T > / < V > for He + using both the DVR and the 
variational method with grid spacing A = 0.4oo. t xy, yz and 

xz components * 2 z 2 - x 2 - y 2 and x 2 — y 2 components 


-•-Iso -»-2po 2pii ^-2so 
_ g r u r u _g 



FIG. 1: The relative error in the n = 1,2 energies of 
vs. internuclear distance, with A =0.25, 0.5, and 0.8ao, from 
R = 0 to 5ao- 


the Coulomb matrix elements are evaluated exactly). We 
have not been able to test the variational method for a 
two-electron problem with the methods available to us, 
due to the prohibitive cost of computing and storing ma¬ 
trix elements. Here we compare the results for the hy¬ 
drogen atom, which tests the one-electron operator. 

In Table U we show the n = 1, 2, 3 eigenvalues of He + 
and the ratios ( T) /( V) for both the DVR and variational 
methods, using a grid spacing of A = 0.4ao- The present 
DVR method clearly gives better energies, in one case 
(2p) by nearly two orders of magnitude. Results for the 
virial theorem are even more decisive, up to three or¬ 
ders of magnitude. We speculate that the reason for this 
favorable performance is that the relationship between 
the Coulomb potential and the kinetic energy operator is 
maintained in matrix form. 



FIG. 2: The relative error in the ground state energy of He 
vs. the number of grid points per spatial dimension with a 
fixed box size of 3ao, plotted on a log-log scale. 


Relative Error in H 2 Ground State Energy 



Internuclear Distance (bohr) 


FIG. 3: The relative error in the ground state energy of H 2 
vs. internuclear distance, with A =0.5 and 0.8ao, from R = 0 
to 5ao, defined relative to the benchmark results of Ref. [4l| . 


B. H+ 

In Figure [T] we show the relative error in the n = 1,2 
energies of H^ for different grid resolutions, compared 
with exact results obtained in prolate spheroidal coordi¬ 
nates as in Ref. [|fj. The errors are on the order of a 
millihartree for A = 0.5 and 0.8ao and are 1-2 orders of 
magnitude better with A = 0.25ao- It appears that the 
2pn u state is generally the most accurate. The errors are 
also relatively constant with respect to the internuclear 
distance. 


C. Two-electron results 

In Figure O we plot the relative error in the ground 
state of Helium for multiple grid resolutions with a fixed 
box size of 3ao, which is sufficient to eliminate truncation 
error. The figure demonstrates a roughly quadratic con¬ 
vergence rate of the ground state energy of Helium with 
respect to the grid resolution; the error is proportional 
to N p where p ~ —2.13. 

In Figure [3] we show the relative error in the ground 
state energy of H 2 for different grid resolutions. As ex¬ 
pected, the results with A = 0.5ao are more accurate 
than with A = 0.8ao. However, both resolutions have 
error roughly on the order of 10 , with the errors being 

slightly larger for R close to zero. Comparing the ground 
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T=1000 
F.T. (arb.units) 



4 6 8 10 12 14 

Energy (eV) 



<3 5 


E U) 
(/) 
CD O 



Energy (eV) 


FIG. 4: Bound-state absorption spectrum for methane calcu¬ 
lated one pulse and equal wall clock computation time. Top 
to bottom: Grid spacing 0.39495414ao, 63 points on a side; 
grid spacing 0.39495414ao, 31 points on a side; grid spacing 
0.19747707ao, 63 points on a side. Time is in atomic units; 
divide by 41.34 for femtoseconds. 


state errors of H 2 with those of H , we see that the H 2 
calculations are slightly less accurate, due to the error 
introduced by the two-electron operator. 


D. Methane photoexcitation, time-dependent 
calculation 


We calculate time-dependent nine-orbital full configu¬ 
ration interaction electronic wave functions for methane 
using the method described in Ref. [42j]. We use grid 
spacings A of 0.39495414 and 0.19747707 bohr. These 


spacings permit a bond length of of 1.086A, which is the 
equilibrium bond length at the highest level of theory 
in Ref. [43|. We position hydrogen nuclei at ( x,y,z ) = 
{(A, A, A), (A,-A,-A), (-A, A,-A), (-A , -A, A)} 
where A = 3A or A — 6A, respectively, for the two grid 
resolutions. The absorption spectra are calculated in a 
straightforward manner, as in Ref. [44j], and shown in 
Fig. HI 

We perform three nine-orbital full-configuration- 
interaction calculations - 15876 Slater determinants, 
5292 spin adapted singlet configurations - with differ¬ 
ent grid bases but otherwise identical, and approximately 
equal total wall clock computation time. Due to the 
different rates of calculation, different final times are 
reached: 

• Resolution A= 0.39495414ao, 63 points on a side: 
tfinal ~7.5fs (300 atomic units), or 

• Resolution A= 0.39495414ao, 31 points on a side: 
tfinal ~24.2fs (1000 atomic units), or 

• Resolution A= 0.19747707ao, 63 points on a side: 
tfinal «121fs (5000 atomic units) 

The difference in speed between the two 63-point calcu¬ 
lations is due to the different behavior of the method for 
the different resolution grids. The difference between the 
two calculations with coarse resolution is due to the size 
of the basis. The basis size affects the Fourier transform 
time (0(n 3 log(n))) and the kinetic energy time (0(n 4 ), 
but with a smaller coefficient). 

Excitation of the low-lying valence states leads to pho¬ 
todissociation; photodissociation of methane has been 
studied by several authors |43|, l45l - l47| . The vertical 
excitation energy of the lowest state was calculated in 
Ref. [43| to be 10.60eV and this appears to be the most 
reliable value in the literature. We have used the same 
bond length as in that work. 

One can see that, of the three panels shown in Fig. HI 
the first one shows a peak that is just above lleV, 
whereas the peak in the other panels occurs higher, be¬ 
tween 12 and 13eV. The difference is due to the trunca¬ 
tion error, the error caused by insufficient spatial extent 
of the grid; the grid used for the top panel is about twice 
as wide as that used for the other two. The difference be¬ 
tween the middle and bottom panels in this figure is due 
to resolution. The grids have approximately the same 
extent in these figures, but the bottom figure uses a grid 
spacing A fti0.2ao, whereas the middle uses A «0.4ao. 
The difference in the t = 300 curves on the middle and 
bottom panels is minor. So we see that a box size of 
24ao and a grid spacing of 0.4ao are probably sufficient 
for calculations of excited state physics on methane. 

The fact that we must use a large basis set to repre¬ 
sent electronic wave functions with large spatial extent 
is certainly a problem with the representation, and we 
comment on the possibility of stretching the grid in the 
asymptotic region later, in the conclusion. We would 
like to avoid truncation error, what comes from a grid 























of insufficient spatial extent, in the analysis of this basis 
set method, and focus on resolution error. In subsec¬ 
tions 1III FI and IIIIGl we perform apples-to-apples com¬ 
parisons of the sine DVR to Gaussian basis sets, in which 
we eliminate truncation error, which does not apply to 
Gaussian basis sets, from consideration. The method we 
use for this purpose is described below. 


E. Method for extrapolation to infinite grid size 

Below we report transition energies and ionization po¬ 
tentials for polyatomic molecules. The desire is to present 
the results without error due to the use of a grid with fi¬ 
nite extent. The “truncation error,” the error due to 
having an insufficient number of points on a side n, is 
uninteresting and should be eliminated. 

The resolution error, in contrast, is of prime im¬ 
portance. We have conjectured that this is an ideal 
smoothed Coulomb representation on a Cartesian grid. 
So we report numbers that are functions of resolution A, 
but not of points on a side n. However, they have error 
bars due to the method used to extrapolate them, as a 
function of points on a side n, to n = oo. We report ex¬ 
citation energies as a function of resolution A, with error 
bars due to truncation error in our finite basis calcula¬ 
tions. The method we use to extrapolate the energies is 
ad hoc and is as follows. 

We choose a function f(n;P) that is monotonic as a 
function of n and that approaches a limit as n —> oo 
but that is otherwise arbitrary, and that is a function not 
only of n but also of certain number Np of parameters 
Pi, i = 1 ...Np. Each eigenvalue E a (n ), a = 1... 12 is fit 
separately as a function of n to the function / by varying 
the parameters of /. 

The uncertainty in the extrapolated energy eigenvalue 
will be affected by the choice of /. We regard / as un¬ 
known and seek to find one that provides acceptable pre¬ 
cision in the reported extrapolated eigenvalue. Presently 
we have tried functions of the form Pi + P 2 rfi e~ nP3 and 
find that 


/(n; Pi, P 2 , P 3 ) = Pi + P 2 n 2 e~ nP3 (36) 


gives a consistently superior fit to the present data, when 
compared to the other choices we tried with Q ^ 2, so 
we chose this function, with an n 2 factor in the exponen¬ 
tial term, for /. We perform a least squares regression, 
choosing Nc values of n with which to perform calcula¬ 
tions and minimizing 


7a=1...12,i=l„ JVp 


dP n 


No 2 

El — f (?lj 1 Pa ) 


3 =1 


= 0 


(37) 

The predicted asymptote is the first parameter, Pi from 
equation 1361 the constant term, 


(38) 


The variance in the predicted E a ( 00 ) will be denoted 
<t 2 . There is systematic error in the prediction due to the 
lack of knowledge about the exact form of the unknown 
function /. There is statistical error due to imperfect 
convergence of the calculated eigenvalues E a (n). So we 
estimate the variance as 

al = {aTf + K at f (39) 

with the systematic error defined as the asymptotic stan¬ 
dard error of the parameter Pi. 

The statistical error for each computed eigenvalue 
E a (n) is that caused by imperfect convergence of the 
MCTDHF relaxation procedure. We have a primitive 
implementation but choose a stringent convergence cri¬ 
terion. The change in energy between the penultimate 
and final iterations, which we will denote A E a (n), is 
generally less than one microhartree, and this number is 
recorded for each eigenvalue and used to estimate a^ at . 
We performed several small runs with an error criterion 
even more stringent. We estimate that the change in 
energy between the penultimate and final iterations is 
significantly more than 100 times the error in the final 
eigenvalue, and therefore we conservatively estimate the 
statistical error for each point separately as 

<j s a tat {n ) = 100 x A E a (n) . (40) 

Given that these individual statistical errors may be cor¬ 
related, we define the statistical error of the overall fit as 
the average of them, 


, n c 

^ = FEOi) (41) 

Nc U 

In summary, we conservatively define the variance in 
the fitted asymptote, the variance in fitted value of the 
transition energy in the limit of infinite basis size, as 

a 2 = {aTf + ^ £ AE ^) j ( 42 ) 

Furthermore, we perform two calculations with differ¬ 
ent values of the parameters Ubi g and n sma ii in order to 
check the error due to the approximations made in our 
calculation of the Coulomb matrix elements. The signifi¬ 
cant figures reported in sections IIII FI and IIII Gl agree for 
the two choices ( nbig,n sma u) = (248,31) and (195,39). 


F. Methane excitation energies 

We calculate excitation energies of methane using the 
same nine-orbital full-configuration-interaction represen¬ 
tation used for the time-dependent calculations above, 
using two The calculation we perform is called state- 
averaged multiconfiguration self-consistent field (MC- 
SCF) and consists of minimizing the average energy of 


E a { 00 ) = Pi 
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A = 

0.39495414a o 

0.19747707 


DZ-s 

DZ 

TZ-s 

TZ 

QZ-s 

QZ 

T 

9.666150(l)eV 9.5353(5) 

T 

9.5845 

9.5828 

9.5381 

9.5359 

9.5197 

9.5160 

E 

10.718654(3) 

10.6309(2) 

T 

10.8820 

10.8815 

10.7792 

10.7774 

10.7462 

10.7419 

T 

10.772476(3) 

10.6826(2) 

E 

10.9799 

10.9795 

10.8225 

10.8173 

10.7585 

10.7517 

T 

10.773745(3) 

10.6850(2) 

T 

11.0339 

11.0335 

10.8792 

10.8740 

10.8157 

10.8089 


TABLE II: Transition energies E a ( oo) for methane calculated 
with 12-state-averaged MCSCF using the sine DVR basis, for 
two grid resolutions, calculated by extrapolating to infinite 
basis size using the method of Sec. IIII El 


TABLE III: 12-state-averaged MCSCF energies calculated 
with Gaussian basis sets, aug-cc-pvdz, aug-cc-pvtz, and aug- 
cc-pvqz, either the full Cartesian basis (no extension) or con¬ 
tracting them spherically (extension -s above). 


the first twelve electronic states of methane with respect 
to variations both of the coefficients of the sine DVR basis 
functions comprising the nine orbitals, and of the coeffi¬ 
cients of the spin-adapted linear combinations of Slater 
determinants. 

These energies are calculated as a function of grid res¬ 
olution, independent of box size (points on a side n), but 
with error bars that are due to finite box size calcula¬ 
tions, using the method described in the subsection im¬ 
mediately above, and reported in Table El Nine or eight 
calculations are used for the extrapolation to infinite ba¬ 
sis size, respectively: for A=0.39495414ao, n =105, 115, 
. . . 185; for A=0.19747707, n =135, 145, 155, . . . 215. 

For comparison, we perform the same state-averaged 
MCSCF calculations using Gaussian basis sets, using the 
Columbus suite of codes for quantum chemistry [H]. We 
use three basis sets, aug-cc-pvdz, aug-cc-pvtz, and aug- 
cc-pvqz (4)1, using either the full set of Cartesian basis 
functions or contracting them to make spherical harmon¬ 
ics. These results are reported in Table uni 

There are only two columns in Table El for only 
one molecule; any conclusions about the method at this 
stage must be considered preliminary. The columns in 
Table El the results with the current sine basis, dif¬ 
fer consistently by about O.leV. The double-zeta and 
triple-zeta columns in Table Ell obtained with standard 
Gaussian basis set methods, have a range of differences, 
from 0.05 to O.lGeV. Therefore, it appears that rougly 
double-zeta accuracy is obtained with a grid spacing 
A=0.39495414ao, and roughly triple-zeta accuracy is ob¬ 
tained with A=0.19747707. We perform a more quanti¬ 
tative analysis of the performance of the representation 
as a function of grid resolution A in the next section, on 
cubane. 


G. Cubane (CsHs) ionization potential 

As presently described, without elaboration, this rep¬ 
resentation for electronic wave functions of molecules us¬ 
ing the sine discrete variable representation (DVR) re¬ 
quires that nuclei be placed on the Cartesian grid points 
and as such, has limited applicability. In the conclu¬ 
sion we speculate about elaborations to the method that 
would allow it to calculate a molecule in an arbitrary 
internuclear geometry. 


For the moment, the cubane molecule provides a good 
test of the method due to its cubic geometry. Not only 
is it cubic, but the C-C and H-H distances are approxi¬ 
mately in the ratio 9:5 or 1.8. The theoretical equilibrium 
geometry calculated at the coupled cluster with single 
and double excitations (CCSD) using the cc-pVDZ Dun¬ 
ning basis set [49j, as tabulated by NIST 0], has the car¬ 
bons at x,y,z = ±0.7893 Angstrom and the hydrogens at 
1.4248 Angstrom, a ratio of 1.805. So we take the geomet¬ 
ric average of these distances, and multiply and divide by 
the square root of 1.8, to arrive at our geometry, with the 
carbons at x,y,z = ±1.493691ao (approximately 0.7904 
Angstrom) and the hydrogens at x,y,z = ±2.6886438ao 
(approximately 1.4228 Angstrom). We use three grid res¬ 
olutions, A = 0.5974764, 0.2987382, and 0.1493691a o . 

In Table E0 we present results showing the first two 
ionization potentials of cubane in the Hartree-Fock ap¬ 
proximation, calculated as in Sec. IIII El extrapolated to 
infinite basis size, for the three grid spacings A. These 
infinite-basis results are then extrapolated to A = 0, and 
that result is shown in the fourth row of the table. The 
method that we use for this final extrapolation is de¬ 
scribed later in this section. 

Two potentials are reported, both those corresponding 
to the difference between the fully converged Hartree- 
Fock energies of the neutral and cation, labeled “I.P.” in 
Table ITVl and those corresponding to the Koopman’s ion¬ 
ization potential, labeled “K.I.P.,” corresponding to the 
neutral Hartree-Fock highest occupied molecular orbital 
energies, calculated as the difference between the neu¬ 
tral Hartree-Fock energy and the cation energy obtained 
through diagonalization using the neutral Hartree-Fock 
orbitals. The precision obtained in the latter is much 
lower than the former due to the primitive Hartree-Fock 
implementation we use. Five or six points are used for 
the extrapolation to infinite basis size; for resolution A = 
0.5974764a 0 , n =64, 72, 80, 90, 108; for A = 0.2987382, 
n =81, 91, 99, 105, 117; and for A = 0.1493691a o , 
n =185, 195, 205, 215, 225, and 235. 

However, the precision in the results for cubane in Ta¬ 
ble IIVI does not come from the extrapolation. For the 
ionization potentials (I.P.) the precision comes from dis¬ 
agreement between the two calculations for the different 
choices for n^ig and n sma u ; for the Koopman’s ionization 
potentials (K.I.P.) the precision comes from nonconver¬ 
gence of the primitive Hartree-Fock procedure, and our 
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conservative choice for the definition of statistical error 
based upon it. 

The ionization potentials of cubane have been previ¬ 
ously calculated in Refs. [5lU53j . The lowest T 2g and 
T 2u Koopmans’ ionization potentials, exactly analogous 
to those reported here, were calculated to be 10.39 and 
10.58eV at the double zeta with polarization level of the¬ 
ory, in Ref. [5^]. In a different basis, the K.I.P.s were 
10.42 and 10.59eV, and the delta-SCF result, closely com¬ 
parable to the I.P. reported here, was 9.74eV for both 
the T 2g and T 2u states. In Ref. j52[, the Koopmans’ I.P.s 
were calculated as 10.40 and 10.62eV, respectively, and 
the I.P.s were calculated to be 9.38 and 9.73eV at a higher 
level of theory with more correlation. 

By comparing these numbers from the literature to 
those in Table IIV1 it seems that the present represen¬ 
tation will be able to produce qualitatively accurate re¬ 
sults on polyatomic molecules using a grid spacing of 
approximately 0.3ao, the medium resolution in the ta¬ 
ble. At this medium resolution, it seems that double-zeta 
quality transition energies are obtained; stepping up to 
the finest resolution produces improvements of less than 
O.leV for the lower cation ionization potential, and im¬ 
provements slightly greater for the higher I.P. This ac¬ 
curacy of O.leV is unsatisfactory for many applications 
involving ground-state Born-Oppenheimer dynamics; the 
standard called for there, “chemical accuracy,” is one 
kilocalorie per mole [35j , which is approximately 43meV. 
Examining the lowest-resolution results, one can see that 
errors introduced going from the resolution of approxi¬ 
mately A = 0.3ao to 0.6ao are a substantial fraction of 
an electron volt. The accuracy at 0.6ao is probably unsat¬ 
isfactory for almost all applications, but A = 0.3ao seems 
sufficient for qualitative studies of excited state potential 
energy curves and time-dependent electron dynamics of 
polyatomic molecules. 

We have attempted to quantify the performance more 
accurately by extrapolating the results in the first three 
rows of Table EY] to A = 0. The power law for the er¬ 
ror that was observed for the one electron results and re¬ 
ported above - a power law A 213 for the error does not 
fit the results on cubane in Table HVl Unfortunately, the 
exponent in the power law for these cubane results is sig¬ 
nificantly lower. In order to obtain error bars on the pre¬ 
dicted extrapolation, we fit the four columns of Table IIVI 
to the same power law. In other words we consider the 
columns in the table to be labeled Ei( A), and with these 
twelve points fit the nine parameters {ei...e 4 , & 1 ...& 4 , Q} 
in the functional form 


E i {A) = e i + b i AQ (43) 

using this fit we obtain the power law exponent Q = 
1.205 ± 0.18. The error bars in the final row of Table HVl 
showing this extrapolation, are almost entirely due to the 
error of this fit, and not to the error of the points used 
in the fit. 

Using the values of bi from this fit, we obtain the 
value at which the accuracy of the computed ioniza- 


Resolution 

2 T 2g 

I.P. K.I.P. 

2m 

J- 2u 

I.P. K.I.P. 

0.5974764 

9.84027(1) 

10.68006(1) 

10.24536(2) 

11.16661(2) 

0.2987382 

9.69775(2) 

10.54452(50) 

9.84952(2) 

10.79933(40) 

0.1493691 

9.63852(10) 

10.47702(10) 

9.69365(4) 

10.64553(4) 

0 

9.591(2) 

10.433(5) 

9.560(9) 

10.523(3) 


TABLE IV: Ionization potentials of cubane, in electronvolts, 
as described in the text. 


tion potentials for cubane is one kilocalorie per mole 
or 43meV, also known as “chemical accuracy” [35|. By 
solving 0.043 = feiA 1 205 , in electronvolts, for A, given 
the fitted bi , we obtain A = 0.139ao from both the I.P. 
and the K.I.P. of the lower ( 2 T 2g ) state, and approxi¬ 
mately A = 0.06cio for both I.P. and K.I.P. for the up¬ 
per ( 2 T 2 „) state. The average of these values is about 
0.1 bohr. Given the flexibility of the representation, it 
seems reasonable to expect that chemical accuracy will 
be obtained generally, for other molecules as well, at this 
resolution. In the conclusion, we mention improvements 
to the method that would account for the truncation of 
the basis in momentum space and that would hopefully 
yield chemical accuracy with an even lower resolution. 


IV. CONCLUSION 


We have demonstrated an efficient real-space basis set 
representation for electronic structure using sine basis 
functions, a generalization of the method of Ref [17] to 
Cartesian coordinates. This and that method make use 
of a resolution-of-the-identity approximation to arrive at 
diagonal expressions for the one- and two-electron matrix 
elements. The singular Coulomb potential is discarded 
and the one- and two-electron matrix elements are ob¬ 
tained instead from the kinetic energy matrix elements 
by requiring that the relationship between the Coulomb 
potential and Laplace operator - that the former is the 
Green’s function of the latter - be maintained in their 
numerical matrix representations. The normally forth- 
rank tensor of two-electron matrix elements is rendered 
first-rank, and may be stored in memory for even the 
largest problems. The energies and virial theorem ratios 
calculated are far superior to those obtained with the 
variational method using the same sine basis. 

We note that this DVR representation bears similar¬ 
ity to that of Ref. [54j, a three-dimensional treatment for 
atoms in spherical coordinates. The three dimensional 
representations, the present one and that of Ref. ]54j , 
as opposed to the treatment for spherical coordinates in 
Ref. [ 13 ] in which only the radial coordinate is discretized, 
permit the maximum degree of parallel computer scala¬ 
bility. We also note that a similar ansatz involving the 
kinetic energy operator has been applied to Gaussian ba¬ 
sis functions in Ref. (55j . 

The representation described in this paper may pro- 
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vide a foundation for an efficient treatment of electronic 
structure that would compete with Gaussian basis set 
methods in applications for which chemical accuracy [35j 
is required. With this goal in mind, several easy-to- 
implement elaborations that would improve its perfor¬ 
mance are conceivable. For instance, it is desirable to 
have grids with different resolutions for different elec¬ 
trons in the Slater determinant basis, such that different 
orbitals with different spatial extents can be described 
efficiently. Algebra along these lines is presented the Ap¬ 
pendix. Also, effective theory [56} may be used to account 
for the truncation in momentum space and improve the 
convergence of the results with respect to resolution A. 

However, the most important improvement is to per¬ 
mit small grid distortions. Presently, the method requires 
that nuclei only be placed on the Cartesian grid points, 
which is its most major limitation. If the grid could be 
distorted slightly, but arbitrarily, then arbitrary internu- 
clear geometries could be calculated simply by distort¬ 
ing the grid such that the grid points and nuclei coin¬ 
cide. Furthermore, if these distortions could be made 
complex-valued, then the representation would be capa¬ 
ble of calculating ionization using the method of complex 
coordinate scaling [571461} . 

Implementing complex-valued grid distortions is there¬ 
fore the next step in the development of this real- 
space representation for electronic structure. Includ¬ 
ing grid distortions in the method will permit arbi¬ 
trary fixed-nuclei geometries and the accurate represen¬ 
tation of ionization. It will also enable calculations of 
fully nonadiabatic electronic and nuclear dynamics of 
polyatomic molecules subject to intense, ultrafast laser 
light, with th e op en-source implementation published in 
Refs. [H, HI, [63f. Because of its efficiency and uniform 
resolution, this sine discrete variable representation for 
electronic structure is best suited to highly correlated, 
highly excited dynamics of electrons in molecules, not 
ground state electronic structure. With the method of 
Domcke and coworkers [63466} . we are using it to calcu¬ 
late phase matched signals for wave mixing experiments 
on polyatomic molecules, and we look forward to this and 
other applications in the future. 
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Appendix A: Derivation for different electron one 
and electron two bases 

It is wasteful to define orbitals all in the same basis 
extending over the entire molecule. Many electrons, no¬ 
tably core electrons, will be localized. To account for this, 
it is imperative to define orbitals on different grids with 
different spatial extent and resolution. Such a treatment 
then calls for Slater determinants belonging to different 
classes containing different numbers of electrons occupy¬ 
ing orbitals belonging to different grids. However, since 
the one electron bases are not combined there is no prob¬ 
lem of linear dependence nor any significant issue related 
to orthogonality. 

One could, for example, interpolate the density on the 
sparser grid onto the finer grid, then use T^ 1 for the finer 
grid to evaluate the integral. However, it is interesting 
to try to adapt the derivation directly to the case of two 
different bases. 

The derivation with two different bases for electrons 
one (ij) and two ( kl ) follows. We define 

y kl (n) = (d 3 r 2 Xk\^)x?\r 2 ) 1 „ (Al) 

J \ri~r 2 \ 

Applying the Laplacian to both sides of equation Eq. Ell) 
and approximating y kl (fi) as a linear combination of the 
basis functions in f 2l 

y k \r i)« E ynX { n\n) , (A2) 

ALL n 

applying the Laplacian with respect to ri to this expan¬ 
sion, multiplying through by Xm (fi), integrating over ri 
and equating the results leads to 

E Vn T mn = 2t t f d 3 f x k ] {r)xf ] (f)%m (r) (A3) 

n 

where T mn = -\{Xm |V 2 |Xn ) are the kinetic energy 
matrix elements. 

Again using the resolution of the identity to approxi¬ 
mate the density (sum of squares of basis functions), the 
right hand side of Eq. (EU1) is evaluated as 

E y k n T mn = 2 ttA~ 3/2 6 kl S km (A4) 

n 

with S the overlap matrix 

Sim = J d 3 r x[ 2) ( r)Xm (r) (A5) 

giving 

y k n = 2 ttA ~ 3 S kl E SmniT- 1 )^. (A6) 
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Using these coefficients and inserting Eas. (IA2l) into [T4l 
gives the expression for the two-electron matrix elements 

[ij\kl] = J2vn [ (A7) 

n 

Once again, the integral on the RHS is evaluated by a 
resolution of the identity, this time in the n density, to 


obtain the expression 

[ij\kl] = 27r(A (1) A (2) )- 3 /%4; S in S mk T-'. (A8) 

mn 

Because the electron one and electron two grids are not 
commensurate, more than one column (equivalently, with 
different indexing, one row) of T~ l will have to be stored. 
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